/*
 * lognumber.cpp
 *
 *  Created on: Mar 9, 2011
 *      Author: tqlong
 */

#include <cmath>
#include <boost/math/special_functions/log1p.hpp>
#include <limits>
#include <algorithm>
#include <cassert>
#include "lognumber.h"

double lognum::DOUBLE_INFINITY=std::numeric_limits<double>::infinity();
double lognum::TOLERANCE=1e-6;

std::ostream& operator<<(std::ostream& s, const lognum& x)
{
  s << (x.sgn()<0?"-":"+") << "exp{" << x.log() << "}";
  return s;
}

namespace std {
  lognum log(const lognum& x)
  {
    assert(x.sgn() >= 0);
    return lognum(x.log());
  }
  lognum log1p(const lognum& x)
  {
    assert(x.sgn() >= -1);
    return lognum(::log1p(::exp(x.log())*x.sgn()));
  }
  lognum sqrt(const lognum& x)
  {
    assert(x.sgn() >= 0);
    return lognum(x.log()/2, x.sgn());
  }
  lognum exp(const lognum& x)
  {
    return lognum(::exp(x.log())*x.sgn(), 1.0);
  }
  lognum logistic(const lognum& x)
  {
    return lognum((-log1p(exp(-x))).toDouble(), 1);
  }
}

lognum::lognum(double d, double sign)
  : log_(d+::log(::fabs(sign))), sign_(sign >= 0 ? 1 : -1)
{
}

lognum::lognum(double d)
  : log_(::log(::fabs(d))), sign_(d >= 0 ? 1 : -1)
{
}

lognum& lognum::operator=(double d)
{
  log_ = ::log(::fabs(d));
  sign_ = d > 0 ? 1 : -1;
  return *this;
}

bool lognum::operator <  (const lognum& x) const
{
  if (log_ == -DOUBLE_INFINITY && x.log_ == -DOUBLE_INFINITY) return false;
  if (sign_ < x.sign_) return true;
  else if (sign_ > x.sign_) return false;
  else if (sign_ > 0)
    return log_ < x.log_ - TOLERANCE;
  else
    return x.log_ < log_ - TOLERANCE;
}

bool lognum::operator >  (const lognum& x) const
{
  return x < *this;
}

bool lognum::operator <= (const lognum& x) const
{
  if (log_ == -DOUBLE_INFINITY && x.log_ == -DOUBLE_INFINITY) return true;
  if (sign_ < x.sign_) return true;
  else if (sign_ > x.sign_) return false;
  else if (sign_ > 0)
    return log_ <= x.log_ + TOLERANCE;
  else
    return x.log_ <= log_ + TOLERANCE;
}

bool lognum::operator >= (const lognum& x) const
{
  return x <= *this;
}

bool lognum::operator == (const lognum& x) const
{
  if (log_ == -DOUBLE_INFINITY && x.log_ == -DOUBLE_INFINITY) return true;
  if (sign_ < x.sign_) return false;
  else if (sign_ > x.sign_) return false;
  else
    return log_ <= x.log_ + TOLERANCE && x.log_ <= log_ + TOLERANCE;
}

bool lognum::operator != (const lognum& x) const
{
  return !(x == *this);
}

double lognum::logsum(double loga, double logb)
{
  double min = std::min(loga, logb),
         max = std::max(loga, logb);
  return max + ::log1p(::exp(min-max));
}

lognum& lognum::operator += (const lognum& v)
{
  if (sign_ == v.sign_) {
    sign_ = sign_;
    log_ = logsum(log_, v.log_);
  }
  else {
    if (log_ >= v.log_) {
      sign_ = sign_;
      log_ = log_ + ::log1p(-::exp(v.log_-log_));
    }
    else {
      sign_ = v.sign_;
      log_ = v.log_ + ::log1p(-::exp(log_-v.log_));
    }
  }
  return *this;
}

lognum  operator +  (double d, const lognum& x) { return lognum(d)+x; }
lognum  operator -  (double d, const lognum& x) { return lognum(d)-x; }
lognum  operator *  (double d, const lognum& x) { return lognum(d)*x; }
lognum  operator /  (double d, const lognum& x) { return lognum(d)/x; }
